Image and cell-level quality control

The following section discusses possible quality indicators for data obtained by IMC and other highly multiplexed imaging technologies. Due to the complexity of the data, the quality metric range across the image and single-cell levels and the chapter is sectioned as such.

Read in the data

We will first read in the data processed in previous sections. As compensation on all images required too much memory (>45GB), we are not using the compensated images below. For convenience, we will re-set the channelNames to their biological targets.

path <- "/mnt/Multimodal_Imaging_Daria/_Data_Analysis/_tmp_daria/Image_analysis/20220723_Steinbock_Mesmer_IFbased"
images <- readRDS(file.path(path,"data/images.rds")) 
masks <- readRDS(file.path(path,"data/masks.rds"))
spe <- readRDS(file.path(path,"data/spe.rds"))

channelNames(images) <- rownames(spe)

Segmentation quality control

The first step after image segmentation is to observe its accuracy. Without having ground-truth data readily available, a common approach to segmentation quality control is to overlay segmentation masks on composite images displaying channels that were used for segmentation. The cytomapper Bioconductor package supports exactly this tasks by using the plotPixels function.

Here, we select 3 random images and perform image- and channel-wise normalization (channels are first min-max normalized and scaled to a range of 0-1 before clipping the maximum intensity to 0.2).

set.seed(20220118)
img_ids <- sample(seq_len(length(images)), 3)

# Normalize and clip images
cur_images <- images[img_ids]
cur_images <- normalize(cur_images, separateImages = TRUE)
cur_images <- normalize(cur_images, inputRange = c(0, 0.2))

plotPixels(cur_images,
           mask = masks[img_ids],
           img_id = "image_id",
           missing_colour = "white",
           colour_by = c("CD8a", "CD3", "CD4", "GD2", "DAPI"),
           colour = list(CD8a = c("black", "yellow"),
                         CD3 = c("black", "red"),
                         CD4 = c("black", "green"),
                         GD2 = c("black", "cyan"),
                         DAPI = c("black", "blue")),
           image_title = NULL,
           legend = list(colour_by.title.cex = 0.7,
                         colour_by.labels.cex = 0.7))

To zoom into the image you can right click and select Open Image in New Tab. We can see that nuclei are centered within the segmentation masks and all cell types are correctly segmented A common challenge here is to segment large (e.g. GD2 - in cyan) versus small (e.g. T cells - in red). However, the segmentation approach here appears to correctly segment cells across different sizes.

An additional approach to observe cell segmentation quality and potentially also antibody specificity issues is to visualize single-cell expression in form of a heatmap. Here, we sub-sample the dataset to 2000 cells for visualization purposes and overlay the sample from which the cells were extracted.

library(RColorBrewer)

n <- length(unique(spe$sample))
qual_col_pals = brewer.pal.info[brewer.pal.info$category == 'qual',]
col_vector = unlist(mapply(brewer.pal, qual_col_pals$maxcolors, rownames(qual_col_pals)))

metadata(spe)$color_vectors$sample <- setNames(sample(col_vector, n), 
                unique(spe$sample))
library(dittoSeq)
library(viridis)

cur_cells <- sample(seq_len(ncol(spe)), 2000)

dittoHeatmap(spe[,cur_cells], genes = rownames(spe)[rowData(spe)$use_channel],
             assay = "exprs", cluster_cols = TRUE, scale = "none",
             heatmap.colors = viridis(100), annot.by = "sample",
             annotation_colors = list(sample = metadata(spe)$color_vectors$sample))

dittoHeatmap(spe[,cur_cells], genes = rownames(spe)[rowData(spe)$use_channel],
             assay = "exprs", cluster_cols = TRUE, scaled.to.max="TRUE",
             heatmap.colors = viridis(100), annot.by = "sample",
             annotation_colors = list(sample = metadata(spe)$color_vectors$sample))

We can differentiate between B-cells (CD20) and T cells (CD3). Some of the markers are specifically detected (e.g., Ki67, CD20, GD2) and others are more broadly detected (e.g. CD45, HLA-ABC, MPO).

Image-level quality control

Image-level quality control is often performed using tools that offer a graphical user interface such as QuPath and FIJI. Viewers that were specifically developed for IMC data can be seen here. In this section, we will specificaly focus on quantitative metrics to assess image quality.

It is often of interest to calculate the signal-to-noise ratio (SNR) for individual channels and markers. Here, we define the SNR as:

\[SNR = I_s/I_n\]

where \(I_s\) is the intensity of the signal (mean intensity of pixels with true signal) and \(I_n\) is the intensity of the noise (mean intensity of pixels containing noise). Finding a threshold that separates pixels containing signal and pixels containing noise is not trivial and different approaches can be chosen. Here, we use the otsu thresholding approach to find pixels of the “foreground” (i.e., signal) and “background” (i.e., noise). The SNR is then defined as the mean intensity of foreground pixels divided by the mean intensity of background pixels. We compute this measure as well as mean signal intensity per image as well as the 95% confidence interval. The plot below shows the
average SNR versus the average signal intensity across all images.

library(tidyverse)
library(ggrepel)
library(EBImage)

image_subset <- getImages(images, 28:30)
cur_snr <- lapply(image_subset, function(img){
    mat <- apply(img, 3, function(ch){
        # Otsu threshold
        thres <- otsu(ch, range = c(min(ch), max(ch)))
        # Signal-to-noise ratio
        snr <- mean(ch[ch > thres]) / mean(ch[ch <= thres])
        # Signal intensity
        ps <- mean(ch[ch > thres])
        
        return(c(snr = snr, ps = ps))
    })
    t(mat) %>% as.data.frame() %>% 
        mutate(marker = colnames(mat)) %>% 
        pivot_longer(cols = c(snr, ps))
})

cur_snr <- do.call(rbind, cur_snr)

cur_snr %>% 
    group_by(marker, name) %>%
    summarize(mean = mean(value),
              ci = qnorm(0.975)*sd(value)/sqrt(n())) %>%
    pivot_wider(names_from = name, values_from = c(mean, ci)) %>%
    ggplot() +
#    geom_errorbar(aes(y = log2(mean_snr), xmin = log2(mean_ps - ci_ps), 
#                      xmax = log2(mean_ps + ci_ps))) +
#    geom_errorbar(aes(x = log2(mean_ps), ymin = log2(mean_snr - ci_snr), 
#                      ymax = log2(mean_snr + ci_snr))) +
    geom_point(aes(log2(mean_ps), log2(mean_snr))) +
    geom_label_repel(aes(log2(mean_ps), log2(mean_snr), label = marker)) +
    theme_minimal(base_size = 15) + ylab("Signal-to-noise ratio") +
    xlab("Signal intensity")

We observe ELAVL4, LUM and S100B to have high SNR but low signal intensity meaning that in general these markers are not abundantly expressed. The Iridium intercalator has high signal intensity but low SNR. This might be due to staining differences between individual nuclei where some nuclei are considered as background. We do however observe high SNR and sufficient signal intensity for the majority of markers. Results are not representative as only images [28:30] from sample 15-1320 were chosen due to memory issues.

Another quality indicator is the image area covered by cells (or biological tissue). This metric identifies ROIs where little cells are present, possibly hinting at incorrect selection of the ROI. We can compute the percentage of covered image area using the metadata contained in the SpatialExperiment object:

width_px <- 700
height_px <- 700

colData(spe) %>%
    as.data.frame() %>%
    group_by(sample_id) %>%
    summarize(cell_area = sum(area),
           no_pixels = mean(width_px) * mean(height_px)) %>%
    mutate(covered_area = cell_area / no_pixels) %>%
    ggplot() +
        geom_point(aes(sample_id, covered_area)) + 
        theme_minimal(base_size = 15) +
        ylim(c(0, 1)) + 
        theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 2)) +
        ylab("% covered area") + xlab("")

We observe that the images with low/high cell coverage are indeed from the samples with a low/high cell number. These two images can now be visualized using cytomapper.

# Normalize and clip images
img1="20211230_13-4339_BM_ROI_004" 
img2="20220204_17-0390_BM_ROI_001" 
img3="20211222_03-0313_BM_ROI_005"
img4="20220105_16-3300_BM_ROI_001"
cur_images <- images[c(img1, img4)]
cur_images <- normalize(cur_images, separateImages = TRUE)
cur_images <- normalize(cur_images, inputRange = c(0, 0.2))

plotPixels(cur_images,
           mask = masks[c(img1, img4)],
           img_id = "image_id",
           missing_colour = "white",
           colour_by = c("CD8a", "CD3", "CD4", "GD2", "DAPI"),
           colour = list(CD8a = c("black", "yellow"),
                         CD3 = c("black", "red"),
                         CD4 = c("black", "green"),
                         GD2 = c("black", "blue"),
                         DAPI = c("black", "blue")),
           legend = list(colour_by.title.cex = 0.7,
                         colour_by.labels.cex = 0.7))

Some samples might be oversegmented due to a high number of neutrophils, which have irregular shapes.

Finally, it can be beneficial to visualize the mean marker expression per image to identify images with outlying cell type compositions. This check does not indicate image quality per se but can highlight biological differences. Here, we will use the aggregateAcrossCells function of the scuttle package to compute the mean expression per image. For visualization purposes, we again asinh transform the mean expression values.

library(scuttle)

image_mean <- aggregateAcrossCells(spe, 
                                   ids = spe$sample_id, 
                                   statistics="mean",
                                   use.assay.type = "counts")
assay(image_mean, "exprs") <- asinh(counts(image_mean))

 image_mean%>%dittoHeatmap(genes = rownames(spe)[rowData(spe)$use_channel],
             assay = "exprs", cluster_cols = TRUE, scale = "none",
             heatmap.colors = viridis(100), 
             annot.by = c("sample", "staining_date"),
             annotation_colors = list(sample = metadata(spe)$color_vectors$sample,
                                      staining_batch = metadata(spe)$color_vectors$staining_date),
             show_colnames = FALSE)

We observe extensive biological variation across the images. Some images contain a high fraction of tumor cells (GD2), T cells (CD3) or proliferating cells (Ki67). These differences will be further studied in the following chapters.

Cell-level quality control

In the following paragraphs we will look at different metrics and visualization approaches to assess data quality (as well as biological differences) on the single-cell level.

Related to the signal-to-noise ratio (SNR) calculated above on the pixel-level, a similar measure can be derived on the single-cell level. Here, we will use a two component Gaussian mixture model for each marker to find cells with positive and negative expression. The SNR is defined as:

\[SNR = I_s/I_n\]

where \(I_s\) is the intensity of the signal (mean intensity of cells with positive signal) and \(I_n\) is the intensity of the noise (mean intensity of cells lacking expression). We calculate the SNR and signal intensity by fitting the mixture model across the transformed counts of all cells contained in the SpatialExperiment object.

library(mclust)

set.seed(220224)
mat <- apply(assay(spe, "exprs"), 1, function(x){
    cur_model <- Mclust(x, G = 2)
    mean1 <- mean(x[cur_model$classification == 1])
    mean2 <- mean(x[cur_model$classification == 2])
    
    signal <- ifelse(mean1 > mean2, mean1, mean2)
    noise <- ifelse(mean1 > mean2, mean2, mean1)
    
    return(c(snr = signal/noise, ps = signal))
})
    
cur_snr <- t(mat) %>% as.data.frame() %>% 
        mutate(marker = colnames(mat))

cur_snr %>% ggplot() +
    geom_point(aes(log2(ps), log2(snr))) +
    geom_label_repel(aes(log2(ps), log2(snr), label = marker)) +
    theme_minimal(base_size = 15) + ylab("Signal-to-noise ratio") +
    xlab("Signal intensity")

This analysis is more representative as it considers all cells, instead of cells of just one sample. We see that the nuclear channels have the highest signal intensity, but quite low SNR. CD14, CD3, GD2, CD8a, GZMB, Ki-67 and CD4 are some of the best markers and also have high SNR and intermediate signal intensity here. CD45 and HLA-ABC are expressed on a lot of cells and have a high signal. Some markers make no sense, such as ELAVL4 and PNMT, they are not nearly as good as CD20. HLA-DR should not have such as bad SNR and should also not be expressed on some many cells (as seen in heatmap above). Maybe this is caused by transformation?

Next, we observe the distributions of cell size across the individual images. Differences in cell size distributions can indicate segmentation biases due to differences in cell density or can indicate biological differences due to cell type compositions (tumor cells tend to be larger than immune cells).

colData(spe) %>%
    as.data.frame() %>%
    group_by(sample_id) %>%
    ggplot() +
        geom_boxplot(aes(sample_id, area)) +
        theme_minimal(base_size = 15) + 
        theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 2)) +
        ylab("Cell area") + xlab("")

summary(spe$area)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    5.00   36.00   51.00   55.85   69.00  748.00

The median cell size is 51 pixels with a median major axis length of 9. The largest cell has an area of 748 pixels which relates to a diameter of r round(sqrt(max(spe$area)), digits = 1) pixels assuming a circular shape. Overall, the distribution of cell sizes is similar across images with image Patient4_005 and Patient4_007 showing reduces average cell size. These images contain fewer tumor cells which can explain the smaller average cell size.

We detect very small cells in the dataset and will remove them. The chosen threshold is arbitrary and needs to be adjusted per dataset.

sum(spe$area < 5)
## [1] 0
spe <- spe[,spe$area >= 5]

Another quality indicator can be the number of cells per image divided by the image size. This measure relates to the percentage of image area covered by cells as explained above.

colData(spe) %>%
    as.data.frame() %>%
    group_by(sample_id) %>%
    summarize(cell_count = n(),
           no_pixels = mean(width_px) * mean(height_px)) %>%
    mutate(cells_per_mm2 = cell_count/(no_pixels/1000000)) %>%
    ggplot() +
        geom_point(aes(sample_id, cells_per_mm2)) + 
        theme_minimal(base_size = 15)  + 
        theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 2)) +
        ylab("Cells per mm2") + xlab("")

The number of cells per mm\(^2\) varies across images which also depends on the number of tumor/non-tumor cells, but also the number of neutrophils which are more prone to oversegmentation.

The data presented here were stained on individual slides. Also, samples were pre-processed in different locations. These (and other) technical aspects can induce staining differences between samples or batches of samples. Observing potential staining differences can be crucial to assess data quality. We will use ridgeline visualizations to check differences in staining patterns:

multi_dittoPlot(spe, vars = rownames(spe)[rowData(spe)$use_channel],
               group.by = "year", plots = c("ridgeplot"), 
               assay = "exprs", 
               color.panel = metadata(spe)$color_vectors$year)

We observe variations in the distributions of marker expression across years. These variations arise mainly from sample preparation differences between samples. Some markers seem to be more problematic such as HLA-ABC, CD274, CD24, CD45 and MPO. HLA-ABC and CD274 might be caused by staining artifacts.

##Dimensionality reduction on arcsinh-transformed counts Finally, we will use non-linear dimensionality reduction methods to project cells from a high-dimensional (40) down to a low-dimensional (2) space. For this the scater package provides the runUMAP and runTSNE function. To ensure reproducibility, we will need to set a seed; however different seeds and different parameter settings (e.g. the perplexity) parameter in the runTSNE function need to be tested to avoid interpreting visualization artefacts. For dimensionality reduction, we will use all channels that show biological variation across the dataset. However, marker selection can be performed with different biological questions in mind.

We will parallelize the tasks to speed up the computation. In order to reduce the required memory per CPU, we will delete unused variables.

library(scater)
rm(selected_patients, selected_roi, selecting_staining_dates, staining_batch, staining_date, test, width_px, year)
## Warning in rm(selected_patients, selected_roi, selecting_staining_dates, :
## object 'selected_patients' not found
## Warning in rm(selected_patients, selected_roi, selecting_staining_dates, :
## object 'selected_roi' not found
## Warning in rm(selected_patients, selected_roi, selecting_staining_dates, :
## object 'selecting_staining_dates' not found
## Warning in rm(selected_patients, selected_roi, selecting_staining_dates, :
## object 'staining_batch' not found
## Warning in rm(selected_patients, selected_roi, selecting_staining_dates, :
## object 'staining_date' not found
## Warning in rm(selected_patients, selected_roi, selecting_staining_dates, :
## object 'test' not found
## Warning in rm(selected_patients, selected_roi, selecting_staining_dates, :
## object 'year' not found
rm(images)
gc() #free unused space
##             used   (Mb) gc trigger    (Mb)   max used    (Mb)
## Ncells  10946593  584.7   16738248   894.0   16738248   894.0
## Vcells 522025645 3982.8 8132970981 62049.7 7099076349 54161.7
set.seed(220225)
spe <- runUMAP(spe, subset_row = rowData(spe)$use_channel, exprs_values = "exprs", BPPARAM = MulticoreParam()) 
spe <- runTSNE(spe, subset_row = rowData(spe)$use_channel, exprs_values = "exprs", BPPARAM = MulticoreParam()) 

After dimensionality reduction, the low-dimensional embeddings are stored in the reducedDim slot.

reducedDims(spe)
## List of length 4
## names(4): UMAP TSNE UMAP_raw UMAP_log
head(reducedDim(spe, "UMAP"))
##                                   [,1]       [,2]
## 20211222_03-0313_BM_ROI_005_1 9.459925 -0.8084583
## 20211222_03-0313_BM_ROI_005_2 9.450486 -0.9311345
## 20211222_03-0313_BM_ROI_005_3 9.589887 -0.8861686
## 20211222_03-0313_BM_ROI_005_4 9.542805 -0.9048243
## 20211222_03-0313_BM_ROI_005_5 9.561148 -0.8776727
## 20211222_03-0313_BM_ROI_005_6 9.549628 -0.9011635

Visualization of the low-dimensional embedding facilitates assessment of potential “batch effects”. The dittoDimPlot function allows flexible visualization. It returns ggplot objects which can be further modified.

library(patchwork)
library(viridis)

p1 <- dittoDimPlot(spe, var = "year", reduction.use = "UMAP", size = 0.2, legend.show = FALSE) + 
    scale_color_manual(values = metadata(spe)$color_vectors$year) +
    ggtitle("year on UMAP")
p2 <- dittoDimPlot(spe, var = "year", reduction.use = "TSNE", size = 0.2, legend.size=2) + 
    scale_color_manual(values = metadata(spe)$color_vectors$year) +
    ggtitle("year on TSNE") +
    theme(legend.text=element_text(size=5))

(p1 + p2) 

p3 <- dittoDimPlot(spe, var = "sample", reduction.use = "UMAP", size = 0.2, legend.show = FALSE) +
    scale_color_manual(values = metadata(spe)$color_vectors$sample) +
    ggtitle("sample on UMAP")
p4 <- dittoDimPlot(spe, var = "sample", reduction.use = "TSNE", size = 0.2, legend.size=2) + 
    scale_color_manual(values = metadata(spe)$color_vectors$sample) +
    ggtitle("sample on TSNE") + 
    theme(legend.text=element_text(size=2))

(p3 + p4)

multi_dittoDimPlot(
  spe,
  vars = rownames(spe)[rowData(spe)$use_channel],
  reduction.use="UMAP",
  assay="exprs",
  legend.size=0.2,
  size=0.2,
  ncol=3,
  nrow=13,
  min.color="gray90",
  max.color="red"
)  + 
    theme(legend.text=element_text(size=1))

## NULL

We observe a strong separation of T-cells (CD3+ cells) between the patients. However a seperation of tumor (GD2+, CD56+) and non-tumor (CD45+) cells cannot really be seen. Therefore we will also run the UMAP on the non-transformed (raw) counts.

##Dimensionality reduction on raw counts Now we will try to run the UMAP on the non-transformed raw counts to see the effect of arcsinh-transformation.

set.seed(220225)
spe <- runUMAP(spe, subset_row = rowData(spe)$use_channel, exprs_values = "counts", name="UMAP_raw", BPPARAM = MulticoreParam()) 
multi_dittoDimPlot(
  spe,
  vars = rownames(spe)[rowData(spe)$use_channel],
  reduction.use="UMAP_raw",
  assay="exprs",
  legend.size=0.2,
  size=0.2,
  ncol=3,
  nrow=13,
  min.color="gray90",
  max.color="red"
)  + 
    theme(legend.text=element_text(size=1))

## NULL

We see that the result is much worse!

##Dimensionality reduction on log-transformed counts As a final comparison, we will run the UMAP on log-transformed counts.

set.seed(220225)
spe <- runUMAP(spe, subset_row = rowData(spe)$use_channel, exprs_values = "log", name="UMAP_log", BPPARAM = MulticoreParam()) 
multi_dittoDimPlot(
  spe,
  vars = rownames(spe)[rowData(spe)$use_channel],
  reduction.use="UMAP_log",
  assay="exprs",
  legend.size=0.2,
  size=0.2,
  ncol=3,
  nrow=13,
  min.color="gray90",
  max.color="red"
)  + 
    theme(legend.text=element_text(size=1))

## NULL

As there is no improvement, we will stick with the arcsinh transformation.

Save objects

The modified SpatialExperiment object is saved for further downstream analysis.

saveRDS(spe, file.path(path,"data/spe.rds"))

Session Info

SessionInfo
## R version 4.2.0 (2022-04-22)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.4 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] patchwork_1.1.1             scater_1.24.0              
##  [3] mclust_5.4.10               scuttle_1.6.2              
##  [5] ggrepel_0.9.1               forcats_0.5.1              
##  [7] stringr_1.4.0               dplyr_1.0.9                
##  [9] purrr_0.3.4                 readr_2.1.2                
## [11] tidyr_1.2.0                 tibble_3.1.8               
## [13] tidyverse_1.3.2             viridis_0.6.2              
## [15] viridisLite_0.4.0           dittoSeq_1.8.1             
## [17] ggplot2_3.3.6               RColorBrewer_1.1-3         
## [19] BiocParallel_1.30.3         cytomapper_1.9.1           
## [21] SingleCellExperiment_1.18.0 SummarizedExperiment_1.26.1
## [23] Biobase_2.56.0              GenomicRanges_1.48.0       
## [25] GenomeInfoDb_1.32.3         IRanges_2.30.0             
## [27] S4Vectors_0.34.0            BiocGenerics_0.42.0        
## [29] MatrixGenerics_1.8.1        matrixStats_0.62.0         
## [31] EBImage_4.38.0             
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.4.0              backports_1.4.1          
##   [3] systemfonts_1.0.4         plyr_1.8.7               
##   [5] sp_1.5-0                  shinydashboard_0.7.2     
##   [7] digest_0.6.29             htmltools_0.5.3          
##   [9] magick_2.7.3              tiff_0.1-11              
##  [11] fansi_1.0.3               magrittr_2.0.3           
##  [13] ScaledMatrix_1.4.0        SpatialExperiment_1.6.1  
##  [15] googlesheets4_1.0.0       tzdb_0.3.0               
##  [17] limma_3.52.2              modelr_0.1.8             
##  [19] svgPanZoom_0.3.4          R.utils_2.12.0           
##  [21] svglite_2.1.0             jpeg_0.1-9               
##  [23] colorspace_2.0-3          rvest_1.0.2              
##  [25] haven_2.5.0               xfun_0.32                
##  [27] crayon_1.5.1              RCurl_1.98-1.8           
##  [29] jsonlite_1.8.0            glue_1.6.2               
##  [31] gtable_0.3.0              gargle_1.2.0             
##  [33] nnls_1.4                  zlibbioc_1.42.0          
##  [35] XVector_0.36.0            DelayedArray_0.22.0      
##  [37] BiocSingular_1.12.0       DropletUtils_1.16.0      
##  [39] Rhdf5lib_1.18.2           HDF5Array_1.24.2         
##  [41] abind_1.4-5               scales_1.2.0             
##  [43] pheatmap_1.0.12           DBI_1.1.3                
##  [45] edgeR_3.38.4              Rcpp_1.0.9               
##  [47] xtable_1.8-4              dqrng_0.3.0              
##  [49] rsvd_1.0.5                httr_1.4.3               
##  [51] htmlwidgets_1.5.4         ellipsis_0.3.2           
##  [53] pkgconfig_2.0.3           R.methodsS3_1.8.2        
##  [55] farver_2.1.1              uwot_0.1.11              
##  [57] sass_0.4.2                dbplyr_2.2.1             
##  [59] locfit_1.5-9.6            utf8_1.2.2               
##  [61] labeling_0.4.2            tidyselect_1.1.2         
##  [63] rlang_1.0.4               later_1.3.0              
##  [65] munsell_0.5.0             cellranger_1.1.0         
##  [67] tools_4.2.0               cachem_1.0.6             
##  [69] cli_3.3.0                 generics_0.1.3           
##  [71] broom_1.0.0               ggridges_0.5.3           
##  [73] evaluate_0.16             fastmap_1.1.0            
##  [75] fftwtools_0.9-11          yaml_2.3.5               
##  [77] knitr_1.39                fs_1.5.2                 
##  [79] sparseMatrixStats_1.8.0   mime_0.12                
##  [81] R.oo_1.25.0               xml2_1.3.3               
##  [83] BiocStyle_2.24.0          compiler_4.2.0           
##  [85] rstudioapi_0.13           beeswarm_0.4.0           
##  [87] png_0.1-7                 reprex_2.0.1             
##  [89] bslib_0.4.0               stringi_1.7.8            
##  [91] highr_0.9                 RSpectra_0.16-1          
##  [93] lattice_0.20-45           Matrix_1.4-1             
##  [95] vctrs_0.4.1               pillar_1.8.0             
##  [97] lifecycle_1.0.1           rhdf5filters_1.8.0       
##  [99] BiocManager_1.30.18       jquerylib_0.1.4          
## [101] RcppAnnoy_0.0.19          BiocNeighbors_1.14.0     
## [103] irlba_2.3.5               cowplot_1.1.1            
## [105] bitops_1.0-7              raster_3.5-21            
## [107] httpuv_1.6.5              R6_2.5.1                 
## [109] promises_1.2.0.1          gridExtra_2.3            
## [111] vipor_0.4.5               codetools_0.2-18         
## [113] assertthat_0.2.1          rhdf5_2.40.0             
## [115] rjson_0.2.21              withr_2.5.0              
## [117] GenomeInfoDbData_1.2.8    hms_1.1.1                
## [119] parallel_4.2.0            terra_1.6-7              
## [121] grid_4.2.0                beachmat_2.12.0          
## [123] rmarkdown_2.14            DelayedMatrixStats_1.18.0
## [125] googledrive_2.0.0         Rtsne_0.16               
## [127] lubridate_1.8.0           shiny_1.7.2              
## [129] ggbeeswarm_0.6.0